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Numerical results for the concurrence and bounds on the localizable entanglement 
are obtained for the square lattice spin-1/2 XXZ-model and the transverse field 
Ising-model at low temperatures using quantum Monte Carlo. 
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Oh, While physicists have been dealing with entangled states since the birth of 

\ quantum mechanics, one has only recently focused on quantifying entangle- 

' ment in general states of quantum systems. 

Entanglement of two spin-1/2 spins in pure states is well- understood, however 
there is no unique definition of entanglement in general mixed states. One 
^ ' natural definition of entanglement of two spins in a mixed state is the entan- 

■ glement of formation, It is natural as it counts the minimum number 

of maximally entangled states needed to construct the state using a restricted 
set of rules. Ep is a property of the reduced density matrix for the two spins in 
question and is a convex function, implying that it is not possible to create Ep 
by mixing states classically. A less desirable property of Ep is however that 
it doesn't capture the forms of entanglement as is present in for instance the 
GHZ state. Thus one can argue that Ep, although natural, is in some cases a 
too restrictive measure of entanglement. 

An opposite extreme is the entanglement of assistance, -Ea^], which is quan- 
tified as the maximal two-spin entanglement achievable by performing any 
types of measurements on the rest of the spins. Ea is also a property of the 
reduced density matrix, but is a concave function. It is a useful concept when 
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entanglement concentration is needed, as for instance when building quantum 
repeaters. However to be of practical use the experimenter must be able to 
perform all kinds of measurement types in order to realize the potential of Ea- 

To make this measure more practical the locahzable entanglement, El, was 
recently introduced [3]. It is similar to Ea but only allowing local measure- 
ments on the spins in the environment. El is not a property of the reduced 
density matrix alone and as such can characterize more complicated features 
of the wave-function not captured by two-point correlation functions. From a 
condensed matter viewpoint this form of entanglement may provide a useful 
concept. Traditionally states of many-body systems are characterized by their 
two-point correlation functions, and many important characteristics such as 
phases and phase transitions are obtained from these. However there are im- 
portant cases where exotic topological orders play important roles requiring 
one to go beyond two-point correlation functions. Thus it becomes vital to 
identify quantities that describe these orders. Recently a connection between 
entanglement and topological order was established [4]. It was shown that the 
existence of the string order parameter in the AKLT-model, a topological 
quantity, is a direct consequence of the presence of long-range locahzable en- 
tanglement. While finding the exact magnitude of El obviously require in- 
formation beyond what is obtainable using two-point correlation functions, it 
is quite remarkable and interesting that it is possible to obtain bounds on its 
magnitude using purely two-point correlation functions [3]. 

In this Letter we will obtain numerical bounds on El as well as calculate 
the entanglement of formation for two general and important 2D quantum 
spin-1/2 models: The XXZ-model and the transverse field Ising-model. 

While a number of theoretical studies have been devoted to quantifying entan- 
glement in models of one dimensional quajiitvim magnets[5,6,7,8,9], materials 
behaving as two dimensional quantum magnets have in the last two decades 
been heavily studied because of their connection to high-temperature super- 
conductivity. Having such well characterized materials makes them interesting 
also from an applied quantum information theoretic point of view. While it 
remains an open issue as to what extent quantum spin systems are actually 
useful for building quantum devices, a study of their entanglement properties 
constitutes a prerequisite. 

For a pure state of two spins lip) — a\ ii) + b\ Ti)+c| iT) + c?| TT); the deviation 
from a product state can be simply expressed in terms of the concurrence 
C = 2\ad — bc\. The concurrence is functionally related to the pure state 
entanglement E as defined by the Von Neumann entropy of one subsystem: 
E = — xlogg x—{l—x) log2(l— x), where x = l/2+\/l — In a remarkable 

paper Hill and Wootters[10] showed that the same functional relation holds 
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between Ep and Cp, where Cp is the concurrence generahzed to a mixed state: 
CF = max|o,7Al- V^- VA^-TAlj. (1) 



where Ai > A2 > A3 > A4 are the eigenvalues of pp. p = pij is the reduced 
density matrix of the spins i and j obtained by tracing over all other spins in 
the system 

P..- jE(^>f)'^r®'^f' (2) 



the brackets denote ground state expectation values and a are the Pauli ma- 
trices. Greek (latin) indexes indicate spin (site) labels, pij — a\ ® ajp*af cr] 
is the time-reversed density matrix. 

There are no simple expressions for Ea- However an upper bound on Ea is 
the concurrence of assistance Ca[2] which is defined[ll] as 

Ca = max^pfcC(7rfe) (3) 

k 



where C is the pure state concurrence and pk is the probability of obtaining the 
pure two-spin state tt^ using a particular measurement procedure on all spins 
except spin i and j. The maximum is taken over all possible measurement 
procedures. Ca can be calculated as[ll] 

Ca = \/a^+ \/a^+ '\/a^+ -/aI. (4) 



Ca {Cp) is the smallest (largest) concave (convex) function coinciding with 
the pure state concurrence for pure states[12]. 

Ca is also an upper bound on E^ which is defined as in Eq. (3( with the 
difference that the maximum is taken over local measurement procedures only. 
Note that because El is defined in terms of the pure state concurrence the 
entropy of one subsystem S{pi — TrjPij) which is an upper bound on Ea 
is not necessarily an upper bound on El- A lower bound on El is gotten 
by calculating the largest singular value of the 3x3 matrix of connected 
correlation functions [3]: 

Q?f-KcT^)-{<){^^)- (5) 
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The Hamiltonian for the XXZ-model is 

H^.. = E {- + 4) + M<^3] ' (6) 

where the sum is taken over all nearest neighbor sites on a square lattice 
and units are chosen such that the exchange coupling is unity. Using the 
symmetries of this Hamiltonian the concurrence can be written 

= ^max I 0, 1 {afa^) + (af a[) | - ^(1 + (afaj))^ - (af + ap^j . (7) 

This expression is strictly valid only when there is no spontaneous symmetry 
breaking. However as shown in Ref. [13] the concurrence is unaffected by spon- 
taneous symmetry breaking for the zero-field XXZ-model, thus this expression 
yields the correct result. The concurrence of assistance is 

The (equal-time) correlation functions in Eqs. (5), (7) and (8) are calculated 
using Monte Carlo simulations employing the Stochastic Series Expansion 
technique[14] with directed-loop updates[15] on a L x L square lattice with 
periodic boundary conditions. This technique utilizes a high temperature ex- 
pansion to map the quantum system to a classical vertex model residing on a 
space-time lattice. This mapping has many parallels with the Suzuki-Trotter 
mapping of a quantum d dimensional system to a classical system in o? -|- 1 
dimensions. The efficiency of the method relies on a fast motion thru configu- 
ration space which is achieved by the non-local directed loop updates leading 
to very short autocorrelation times. Correlation functions of observables diag- 
onal in the representation basis can be read off directly from the configuration 
while correlation functions of off-diagonal observables can be read off from 
how the loop propagates in space-time. This is because the tail and the head 
of the loop act as operator insertions in the spin-configuration. This makes 
two-point correlations particularly easy to measure. It is possible to measure 
higher order correlators as well, but at considerably increased computational 
cost. 

Fig. 1 shows Cp plotted as a function of A for different separations between 
sites i and j at a low temperature {(5 — 80). Decreasing the temperature fur- 
ther did not change the results significantly. The nearest neighbor concurrence 
is the largest and peaks at the antiferromagnetic point where it reaches a value 
0.169±0.002 which is in good agreement with the value 0.1694 gotten from the 
ground state energy. The concurrence for bigger distances between the spins 
reaches a maximum close to A = —1 and decreases rapidly with increasing 
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Fig. 1. The concurrence as function of A for different separations iixi'iy) — {jx^jy) 
between sites i and j for a square lattice XXZ-model with 32 x 32 sites. The inverse 
temperature is /? = 80. Inset is a comparison to the ID result (open symbols, 128-site 
chain at (3 = 80) for the linear nearest and next-nearest neighbor Cp where the 2D 
Cp is multiplied by d=2. 

A and separation between i and j. The concurrence between diagonal neigh- 
bors is bigger than that between linear next-nearest neighbors. The general 
behavior resembles that of the concurrence in the XXZ-chain although the 
actual values here are lower. This is expected on grounds of the "monogamy" 
property of entanglement: The more neighbors a spin is entangled with, the 
less entanglement per pair. However even when correcting for the increased 
number of neighbors Cp is lower for the square lattice than for the chain. This 
can be seen from the inset in Fig. 1. 

For A < — 1, the ground state is a product state of spins aligned in the same 
direction, thus Ca = 0, and all forms of entanglement are zero. For |A| < 1 
the expectation value (cr^) = at any finite temperatures, thus Ca = 1. A 
lower bound on El is gotten from Qf? which is the maximum correlation 
function in the regime A < 1. A plot of Qf? as a function of A for different 
site separations is shown in Fig. 2. For big separations Qfj^ behaves as a power 
law at low temperatures and so decay slowly implying the existence of 
over long distances. The lower bound is largest towards the ferromagnetic 
point, where also the decay with distance is very slow. Thus one expect to 
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Fig. 2. Lower bound on El as function of A for different separations between sites 
i and j for a square lattice XXZ-model with 32 x 32 sites. The inverse temperature 
is (3 = 80. 

see the largest El there. This is similar to the ID situation where the lower 
bound decays as \i — j|-arccos(-A)/7r ^ A-dependent prefactor[16], thus 
also decaying slowest at the ferromagnetic point. 

For A > 1 the spin symmetry is broken by an infinitesimal staggered field and 
Ca deviates from unity. It is interesting to ask what Ca is for infinite separation 
between i and j in the broken symmetry phase. Using {(jfcTj) — (af) (a|) in the 
limit — > oo one finds CA{\i—j\ — ^ oo) = I— ml , where = | (a^) \ is the 
staggered magnetization which at the Heisenberg point is ms(A = 1) = 0.614 
and increases towards ms(A — > oo) = 1. It is interesting to note that true 
long-ranged El is not ruled out by this upper bound at the Heisenberg point 
(A = 1) or at any finite A, in contrast to the Ising case (A oo) where 
Ca = 0. The lower bound on El decreases however rapidly to zero as the 
separation \i — j\ is increased. This is a consequence of the finite excitation 
gap in the spectrum for A > 1. 

The Hamiltonian for the transverse field Ising model is 



{hi) ' 



(9) 
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Fig. 3. Nearest neighbor concurrence for the transverse field Ising model on a square 
lattice. The different curves are for different lattice sizes L x L at different temper- 
atures. Ac = 0.32841 is gotten from Ref. [17]. The inset shows a comparison to the 
transverse Ising chain (open triangular symbols). In the inset the 2D concurrence 
is multiplied by a factor d = 2. 

This model exhibits a phase transition at A = Ac. For A <^ Ac the transverse 
magnetic field dominates and the ground state is similar to free spins in a 
magnetic field. For A ^ Ac , the system is like an Ising model without a mag- 
netic field thus having a doubly degenerate ground state. In an experimental 
situation this degeneracy is broken which must be taken into account when 
calculating correlation functions. 



The concurrence can be expressed as 



= -max ( 0, 1 {afa]) + {afa]) \ - ^{1 + {afa^y - {af + a^), 

I - (^f ^1) I - ^(l^^^f^) (10) 



This expression is also valid for A > Ac regardless of symmetry breaking [13]. 

The results for the nearest neighbor concurrence is plotted in Fig. 3. The curve 
resembles largely that seen for the transverse field Ising-chain (seen in the in- 
set) although the values are lower consistent with entanglement monogamy. 
However, in contrast to the one dimensional case where the concurrence peaks 
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Fig. 4. Bounds on the nearest-neighbor El on a 32 x 32 lattice for f3 = 40. Upper 
bounds are indicated with circular symbols while square symbols indicate lower 
bounds. The lines are guides to the eye. The closed symbols are for a system without 
a symmetry breaking parameter, while the open symbols are obtained by adding a 
small longitudinal field {H^ = 0.001). 

well on the high-field side of the transition, Apgg^j^ ~ 0.8Ac[7], the concur- 
rence peaks here much closer to the critical point. The statistical errors makes 
it difficult to assess the exact location of the peak. 

For A < Ac Eq. (8) is valid. However in the symmetry broken phase there are 
no simple explicit expression for Ca so we use Eq. (4) where the eigenvalues are 
gotten from numerically diagonalizing pp taking care to measure also (a^a^) 
correlators in p. 

The bounds on is shown in Fig. (4). For A < Ac the upper and lower bounds 
give relatively tight bounds on E^. Both increase from zero and approach each 
other as A ^ A^. For A > Ac the symmetric ground state approaches the GHZ 
state for A — > oo, thus El ^ 1. In this case the upper and lower bounds almost 
coincides thus giving very tight bounds on El- Results for the experimentally 
more realistic case when the symmetry is broken is also shown. Here the upper 
and lower bounds approaches zero quite rapidly. For other distances between 
the spins, the lower bound drops rapidly to zero away from Ac- The upper 
bound also decreases away from Ac, but not as fast. These results resembles 
those obtained in ID [3]. 

In conclusion the concurrence and bounds on the localizable entanglement was 
calculated numerically for the 2D XXZ-model and the transverse field Ising 
model. While the gross features of the results are similar to results in ID, it is 
interesting to note the slow decay of the localizable entanglement close to the 
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ferromagnetic point in the XXZ-model as well as the location of the maximum 
concurrence in the transverse field Ising model. 
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